cd 'E:\MATLAB scripts\Experimental confidence bands\Rice_Realizations\Rice_Realizations\Realizations'
filelistA = dir('*.txt');
C = [];
SPECTRA = [];
max_time = 25000/1;
for ii = 1:max(size(filelistA));
    if filelistA(ii).isdir ~= true
        fname = filelistA(ii).name;
        fid = fopen(fname);
        out = textscan(fid,'%s %f %f','delimiter',',','headerlines',1);
        fclose(fid);
    end
    datetime = out{1};
    mass = out{2};
    dt1 = datetime(1);%time experiment started
    sec = etime(datevec(datetime),datevec(dt1));
    dt = 1.2;
    maxsec = floor(sec(end));
    t = 1:dt:maxsec;%%%time into run in seconds
    t = t';
    mass_int = interp1(sec,mass,t);
    M = movmedian(mass_int,5);%%%%Has been 5
    Q = (M(2:end) - M(1:end-1))/dt;%Q = mass flux in grams/sec
    Q = [0;Q];
    if ii == 1 | ii == max(size(filelistA));
        disp(mean(Q))
    end
    er = abs((mean(Q) - 0.3577)/0.3577);
    if er <= 0.1;
        t = t(1:max_time);
        Q = Q(1:max_time)*50;
        T = max(t);%Max time
        nt = size(Q,1);%total number of mass measurements
        Qdetrended = Q - mean(Q);%Centering of data
        ne = 1;%number of ergodic realizations from full data set
        nte = floor(nt/ne);
        Qdetrended_e = [];
        for j = 1:ne
            qdetrended_e = Qdetrended(1+j*nte-nte:j*nte);
            Qdetrended_e = [Qdetrended_e qdetrended_e];
        end
        Qdetrended_e = Qdetrended_e';
        %% MTM Spectra
        nw = 4;%number of windows
        nsim = 2;%number of random simulations based on AR(1) model fit, to figure out confidence bands
        t = t(1:nte);
        t = t';
        MTMe = [];
        for j = 1:ne;
            datax = Qdetrended_e(j,:);
            datax = datax';
            [Mspecred,specred,po,fd,fr,tabMC,tabchi,tabtchi,theored,rho]=RedConf(datax,t,dt,nsim,nw);
            MTMe = [MTMe po];
        end
        if ne > 1
            MTMmean = mean(MTMe');%ensemble average data
        end
        if ne == 1
            MTMmean = MTMe';
        end
        period = 1./fd;
        %% Find change pts
        lnperiod = log(period);%convert data over to log space
        lnMTMmean = log(MTMmean);%convert data over to log space

        dp = 0.001;%pick an increment that is fine, relative to data spacing, to interpolate unevenly spaced data onto.
        lnperiod_int = min(lnperiod):dp:max(lnperiod(2:end));%create vector with spacing dp to interpolate onto.
        lnMTMmean_int = interp1(lnperiod(2:end),lnMTMmean(2:end),lnperiod_int);%Do the interpolation
        mean_power = mean(exp(lnMTMmean_int));
        MTMmean = MTMmean./mean_power;

%         FC = findchangepts(lnMTMmean_int,'Statistic','linear','MaxNumChanges',2);
%         C1 = exp(lnperiod_int(FC(1)));
%         C2 = exp(lnperiod_int(FC(2)));
%         C = [C;C1];
        SPECTRA = [SPECTRA MTMmean'];

        %% Plot data
        figure(1)
        loglog(period,MTMmean,'k','linewidth',0.5);
        if ii == 1;
            loglog(period,MTMmean,'g','linewidth',0.5);
        end
        if ii == max(size(filelistA));
            loglog(period,MTMmean,'r','linewidth',0.5);
            figure(5)
            loglog(period,MTMmean,'r','linewidth',0.5);
        end
        ylabel('Spectral Power')
        xlabel('period (sec)')
        axis([2.4 6e3 0.003 8])
        hold on
    end
end
cd 'E:\MATLAB scripts\Experimental confidence bands\Rice_Realizations\Rice_Realizations'
fd = fd(2:end-1);
period = 1./fd;
SPECTRA = SPECTRA(2:end-1,:);
MTMe_sort = sort(SPECTRA,2);
figure(2)
hold on
F = Curl_cb;
F(243,:) = [1 0 1];
nx = size(SPECTRA,2);
for i = size(MTMe_sort,2):-1:1
    s = ceil((i-0.5)/nx*size(F,1));
    s = F(s,:);
    linetrial = MTMe_sort(:,i);
    linetrial = [0.003;linetrial;0.003];
    periodtrial = [6e3;period;1];
    if i == size(MTMe_sort,2)
        basex = [2.4 2.4 6e3 6e3];
        basey = [0.003 8 8 0.003];
        fill(basex,basey,F(end,:),'LineStyle','none');
    end
    if i >= 1
        fill(periodtrial,linetrial,s,'LineStyle','none');
    end
%     if i == 1
%         fill(periodtrial,linetrial,[1 1 1],'LineStyle','none');
%     end
    %fill(period,MTMe_sort(:,i),s,'LineStyle','none');
    set(gca, 'YScale', 'log')
    set(gca, 'XScale', 'log')
    %loglog(period,MTMe_sort(:,i),'color',s,'linewidth',1)
    axis([2 1e4 1e1 1e5])
end
MTMe_mean = median(MTMe_sort,2);
loglog(period,MTMe_mean,'c','linewidth',1);
figure(5)
loglog(period,MTMe_mean,'c','linewidth',1);
figure(2)
set(gca, 'YScale', 'log')
set(gca, 'XScale', 'log')
ylabel('Spectral Power')
xlabel('period (sec)')
colormap(F)
b = colorbar;
b.Label.String = 'Fraction of spectra with power less than';
%% Generate 95% Confidence band
B95 = [];
P_int = 0.5:0.05:1;
Percents = 1:1:size(MTMe_sort,2);
Percents = Percents./size(MTMe_sort,2);
for i = 1:max(size(MTMe_sort))
    Power = MTMe_sort(i,:);
    b95 = interp1(Percents,Power,P_int);
    b95 = b95(10);
    B95 = [B95;b95];
end
loglog(period,B95,'m','linewidth',0.25);
figure(5)
loglog(period,B95,'m','linewidth',0.25);
legend('Single Realization','Average of all Realizations','95th percentile Confidence Band')
figure(2)
axis([2.4 6e3 0.01 7]);%6
disp(size(MTMe_sort,2))
figure(4)
semilogx(period,B95./MTMe_mean,'r')
axis([2.4 6e3 0 4])
ylabel('95th percentile confidence band - mean signal power')
xlabel('period (sec)')
%% Generate data to compare imposed periodicity signals to confidence levels
%cd C:\Work\Projects\Rice_Pile\Spring_2022_Datasets\Datasets_Scripts_Chloe_work\Grid\
cd 'E:\MATLAB scripts\Experimental confidence bands 191122\Grid\Grid'
Gr = load('Grid.dat');
filelistB = dir('*.txt');
Sig_Strength = [];
for jjj = 1:max(size(filelistB));
    if filelistB(jjj).isdir ~= true
        fname = filelistB(jjj).name;
        fid = fopen(fname);
        out = textscan(fid,'%s %f %f','delimiter',',','headerlines',1);
        fclose(fid);
    end
    datetime = out{1};
    mass = out{2};
    dt1 = datetime(1);%time experiment started
    sec = etime(datevec(datetime),datevec(dt1));
    dt = 1.2;
    maxsec = floor(sec(end));
    t = 1:dt:maxsec;%%%time into run in seconds
    t = t';
    mass_int = interp1(sec,mass,t);
    t = t(1:max_time);
    mass_int = mass_int(1:max_time)*50;
    M = movmedian(mass_int,5);%%%%%has been 5
    Q = (M(2:end) - M(1:end-1))/dt;%Q = mass flux in grams/sec
    Q = [0;Q];
    T = max(t);%Max time
    nt = size(M,1);%total number of mass measurements
    Qdetrended = Q - mean(Q);%Centering of data
    ne = 1;%number of ergodic realizations from full data set
    nte = floor(nt/ne);
    Qdetrended_e = [];
    for j = 1:ne
        qdetrended_e = Qdetrended(1+j*nte-nte:j*nte);
        Qdetrended_e = [Qdetrended_e qdetrended_e];
    end
    Qdetrended_e = Qdetrended_e';
    %% MTM Spectra
    nw = 4;%number of windows, has been 4
    nsim = 4;%number of random simulations based on AR(1) model fit, to figure out confidence bands
    t = t(1:nte);
    t = t';
    MTMe = [];
    for j = 1:ne;
        datax = Qdetrended_e(j,:);
        datax = datax';
        [Mspecred,specred,po,fd,fr,tabMC,tabchi,tabtchi,theored,rho]=RedConf(datax,t,dt,nsim,nw);
        MTMe = [MTMe po];
    end
    if ne > 1
        MTMmean = mean(MTMe');%ensemble average data
    end
    if ne == 1
        MTMmean = MTMe';
    end
    period = 1./fd;
    period = period(2:end-1);
    MTMmean = MTMmean(2:end-1);
    %% logbinned intpolated data
    lnperiod = log(period);%convert data over to log space
    lnMTMmean = log(MTMmean);%convert data over to log space
    dp = 0.001;%pick an increment that is fine, relative to data spacing, to interpolate unevenly spaced data onto.
    lnperiod_int = min(lnperiod):dp:max(lnperiod(2:end));%create vector with spacing dp to interpolate onto.
    lnMTMmean_int = interp1(lnperiod(2:end),lnMTMmean(2:end),lnperiod_int);%Do the interpolation
    mean_power = mean(exp(lnMTMmean_int));
    MTMmean = MTMmean./mean_power;
    n = Gr(jjj,1);
    [val,idx]=min(abs(period-n));
    sig_strength = MTMmean(idx)/B95(idx);
    Sig_Strength = [Sig_Strength;sig_strength];
end
cd 'E:\MATLAB scripts\Experimental confidence bands\Rice_Realizations\Rice_Realizations'
figure(3)
x = [0 1];
y = [1 1];
plot(x,y,'k')
hold on
ms = 20;
plot(Gr(41:45,2),Sig_Strength(41:45),'v','MarkerFaceColor', [0.2754 0.5012 0.9659],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms)
plot(Gr(11:15,2),Sig_Strength(11:15),'s','MarkerFaceColor', [0.1517 0.7447 0.9142],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms)
plot(Gr(21:25,2),Sig_Strength(21:25),'p','MarkerFaceColor', [0.1439 0.9268 0.6545],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms)
plot(Gr(31:35,2),Sig_Strength(31:35),'h','MarkerFaceColor', [0.5036 0.9988 0.327],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms)
plot(Gr(6:10,2),Sig_Strength(6:10),'d','MarkerFaceColor', [0.8047 0.9245 0.2046],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms)
plot(Gr(26:30,2),Sig_Strength(26:30),'^','MarkerFaceColor', [0.9954 0.6534 0.1958],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms)
plot(Gr(36:40,2),Sig_Strength(36:40),'>','MarkerFaceColor', [0.9408 0.3557 0.0703],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms)
plot(Gr(1:5,2),Sig_Strength(1:5),'o','MarkerFaceColor', [0.7738 0.1503 0.0115],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms)
plot(Gr(16:20,2),Sig_Strength(16:20),'<','MarkerFaceColor', [0.4796 0.0158 0.0106],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms)
ylabel('Relative spectral signal amplitude')
xlabel('Input signal amplitude/mean input rate')
axis([0 1 0 18])
legend('Threshold','Period = 6sec','Period = 12sec','Period = 24sec','Period = 48sec','Period = 100sec','Period = 250sec','Period = 500sec','Period = 1000sec','Period = 2000sec','Location','northwest')
bb = colorbar;
bb.Label.String = 'Period of input signal';
colormap(turbo)
set(gca,'ColorScale','log')
caxis([2 2000])